#!/usr/bin/python
'''
Evaluate the mapping results.
'''
import sys
import os

# Parameters.
map_file = sys.argv[1]

def read_file(filepath):
	# define data structure.
	data = {}
	
	# Check if zipped.
	zipped = False
	if filepath.count(".gz"):
		zipped = True
	
	# Open file.
	seq = ""
	seqheader = ""
	first = True
	
	if zipped == True:
		fin = gzip.open(filepath,"r")
	else:
		fin = open(filepath,"r")
	for line in fin:
		# Skip comments.	
		if line[0] == "#": continue

		# first time.
		if first == True and line[0] == ">":
			seq = ""
			seqheader = line.strip().replace(">","")
			first = False
			continue		

		# finished seq.
		if line[0] == ">":
			# save data.
			data[seqheader] = seq.strip()

			# Reset variables.
			seq = ""
			seqheader = line.strip().replace(">","")
			continue

		# get seq.
		seq += line.strip()
			
	# Handle last close and return.
	data[seqheader] = seq.strip()
	fin.close()
	return data

def is_int(c):
	try:
		number = int(c)
		return True
	except (ValueError, IndexError):
		return False

def parse_custom_file(file_path):
	# Read in file.
	dat = []
	used = {}
	fin = open(file_path, "rb")
	for line in fin:
		#if len(dat) > 10000: break
		
		# skip comments.
		if line[0] == "#": continue
		
		# tokenize.
		tmp = line.strip().split("\t")
		tname = tmp[0].strip().replace(">","")
		read = tmp[3].strip().replace(">","")
		ident = float(tmp[8].replace("%",""))
		cov = float(tmp[12].replace("%",""))

		
		# Only keep one alignment per read.
		if read in used: 
			# See if beats previous.
			if ident > used[read][0]:
				# Delete previous.
				del dat[ used[read][1] ]
			else:
				# skip this one.
				continue
		used[read] = [ident, len(dat)]
		
		# Save data.
		dat.append( (ident,cov, read) )
		
	fin.close()
	return dat


# Read in maps.
map_dat = parse_custom_file(map_file)

# Setup permeter study.

res = .05
i = 0.5
iness = []
cness = []
while i < 1.0:
	cness.append(i)
	iness.append(i)
	i += res
#iness.reverse()
#cness.reverse()

results = {}
for x in iness:
	results[x] = {}
	for y in cness:
		results[x][y] = [0,0]


# Loop over each entry.
choseni = 0.9
chosenc = 0.7
for entry in map_dat:
	# Calc param.
	i_act = entry[0] / 100.0
	c_act = entry[1] / 100.0
	
	
	
	if i_act < choseni or c_act < chosenc:
		print entry[2]
		
	continue
	
	# Make cumulative counts.
	for i_min in iness:
		for c_min in cness:			
			# Check contamination threshold.
			if i_act >= i_min and c_act >= c_min:
				# Add contamination count.
				results[i_min][c_min][0] += 1
			else:
				# Add HT count.
				results[i_min][c_min][1] += 1
sys.exit()				
# Switch keys back.
iness.reverse()
cness.reverse()

# Print results as matrix.
print "contamination map", len(map_dat)
row = ""
for c in cness:
	row += "\t" + str(c) + "\t"
print row
result = results
ikeys = sorted(result.keys())
for ikey in iness:
	# Write row start.
	row = [str(ikey)]
	
	# Loop over second key.
	for ckey in cness:
		row.append( "%i\t%i" % (result[ikey][ckey][0], result[ikey][ckey][1]) )
	row = '\t'.join(row)
	print row	


